module scene.scene;

import material.material;
import material.texture;
import objects.triangle;
import objects.trianglemesh;
import objects.bezierpatch;
import objects.abstractobject;
import render.render;
import scene.intersection;
import scene.ray;
import scene.box;
import std.math;
import std.random;
import utils.matrix;
import utils.vector2;
import utils.vector3;

 import std.stdio;

//version = OnDemandBIH;
//version = IntersectWald;

//version = OldBIH;
version = NewBIH;

//version = DisableBIH

struct TriData {
    AbstractObject obj;
    int orgIdx;
}

class Scene {
    private {
        Triangle[] mTris;
        TriAccel[] mTriAccel;
        TriData[] mTriData;

        Box mBoundingBox;
        //indices into mTris, used by BIH to sort triangles
        uint[] mTriIndex;
        //node 0 is root node; mBoundingBox is important to get the initial bounds
        BIH_Node[] mBihNodes;
        uint mBihNodeCount;
        //object size; I'm not sure if simply a scalar would be enough, maybe yes
        Vector3f mAvgObjSize;
        bool mNormInterpolate;

        bool mUpdating;
    }

    public Material globalMaterial;

    //const float BIH_GRID_SCALEDOWN = 1.0f/6.0f; //BIH paper, 3.3.2, item 1
    private const uint AXIS_X1 = 0;

    private struct BIH_Node {
        //right-shift to get index value
        final const uint INDEX_SHIFT = 2;
        //bits to encode axis, corresponds to AXIS_X1
        // bit0 bit1 Meaning
        //  0    0    x
        //  1    0    y
        //  0    1    z
        //  1    1    leaf node
        final const uint AXIS_MASK = 3;

        //bits before index value: flags
        //index value itsself: index into bih_nodes
        //(NOTE: it could be "faster" to make it a pointer, but GC wont see it!)
        int index_flags;

        union {
            struct { //leaf node
                int item;  //index of object
                int count; //number of objects, there
            }
            float[2] clip;  //internal node
        }
    }

    //tuneable parameter: max. triangles per BIH leaf node
    //also used at dynamic calculation
    const uint cMaxTrisPerNode = 8;

    public void buildBIH() {
        calculateBoundingBox();

        if (mTriIndex.length < mTris.length) {
            mTriIndex.length = mTris.length;
            foreach (uint index, inout uint n; mTriIndex) {
                n = index;
            }
        }

        //preallocate BIH nodes
        //assume a fully used tree, with one triangle per node
        //(pressimistic assumption)
        uint prealloc = mTris.length*2+1;
        if (mBihNodes.length < prealloc)
            mBihNodes.length = prealloc;

        mBihNodeCount = 1; //include root node

        version (OnDemandBIH) {
            mBihNodes[0].index_flags = -1;
            mBihNodes[0].item = 0;
            mBihNodes[0].count = mTris.length;
        } else {
            recurseBIH(mBoundingBox, 0, 0, mTriIndex.length, uint.max);
        }

        writefln("nodes=%d", mBihNodeCount);
    }

    private void recurseBIH(Box currentBox, uint nodeIdx,
        uint index_start, uint index_end, uint maxdepth)
    {
        //writefln("node: %d", nodeIdx);

        /*debug {
            //slow debug test: all triangles within box?
            for (int i = index_start; i < index_end; i++) {
                auto tri = mTris[mTriIndex[i]];
                //vertex
                for (int t = 0; t < 3; t++) {
                    //axis
                    for (int a = 0; a < 3; a++) {
                        assert(tri.v[t][a] >= currentBox.p[0][a]);
                        assert(tri.v[t][a] <= currentBox.p[1][a]);
                    }
                }
            }
        }*/

        bool nosplit = (maxdepth == 0);
        nosplit |= ((index_end - index_start) <= cMaxTrisPerNode);

        //is it worth to split sth.?
        if (nosplit) {
            //create leaf node
            mBihNodes[nodeIdx].index_flags = -1;// BIH_Node.AXIS_MASK; //leaf node
            mBihNodes[nodeIdx].item = index_start;
            mBihNodes[nodeIdx].count = index_end - index_start;
            return;
        }

        //find longest axis & calc splitting position
        uint axis = currentBox.maxExtent();
        float max1 = currentBox.p[0][axis];     //"left"
        float max2 = currentBox.p[1][axis];     //"right"
        float split = max1 + (max2-max1)/2.0f;  //"middle"

        //find out the two planes max1 and max2; check each triangle
        //these will be the values in BIH_Nodes.clip[2]
        int storeIndex = index_start;
        for (int i = index_start; i < index_end; i++) {
            if (mTris[mTriIndex[i]].isLeft(split,axis,max1,max2)) {
                //swap(i, j)
                auto tmp = mTriIndex[i];
                mTriIndex[i] = mTriIndex[storeIndex];
                mTriIndex[storeIndex] = tmp;
                storeIndex++;
            }
        }
        //split is at [index_start .. storeIndex] and [storeIndex .. index_end]

        debug {
            float old_max1 = max1;
            float old_max2 = max2;
            for (int i = index_start; i < storeIndex; i++) {
                assert((mTris[mTriIndex[i]].isLeft(split,axis,max1,max2)));
            }
            assert(old_max1 == max1 && old_max2 == max2);
            for (int i = storeIndex; i < index_end; i++) {
                assert(!(mTris[mTriIndex[i]].isLeft(split,axis,max1,max2)));
            }
            assert(old_max1 == max1 && old_max2 == max2);
        }

        //adjust the bounding boxes
        //calculate only from split plane candidate, not from actual split planes
        Box b1 = currentBox;
        b1.p[1][axis] = split;//max1;
        Box b2 = currentBox;
        b2.p[0][axis] = split;//max2;

        //debug writefln("%s %s %s %s %s %s %s %s %s", currentBox.p[0], currentBox.p[1], axis, max1, split, max2, storeIndex, index_start, index_end);

        //create node, i.e. fill BIH_Node at nodeIdx with current values

        if (false && (index_start==storeIndex || index_end==storeIndex)) {
            debug writefln("%s %s %s %s %s %s %s %s %s", axis, currentBox.p[0][axis], max1, split, max2, currentBox.p[1][axis], storeIndex, index_start, index_end);
            //two cases:
            //1. ouch, can't split into two nodes
            //2. one box is simply empty
            if (max1 == max2) {
                recurseBIH(currentBox, nodeIdx, index_start, index_end, 0);
            } else {
                //let it create a leaf node (rec.-depth==0)
                //but also adjust bounding box
                Box bx = currentBox;
                if (index_start==storeIndex) {
                    bx.p[0][axis] = max2;
                } else {
                    bx.p[1][axis] = max1;
                }
                recurseBIH(bx, nodeIdx, index_start, index_end, maxdepth-1);
            }
        } else {
            //append two new nodes to the end of the array
            //save index into struct
            uint child_index = mBihNodeCount;
            mBihNodeCount += 2;

            //xxx: this still can happen now, since we also might insert empty nodes!
            assert(mBihNodeCount <= mBihNodes.length, "wrong max BIH nodes assumption");

            mBihNodes[nodeIdx].index_flags = (child_index<<BIH_Node.INDEX_SHIFT) | axis;
            mBihNodes[nodeIdx].clip[0] = max1;
            mBihNodes[nodeIdx].clip[1] = max2;

            //recurseBIH should return the index of the created BIH_Node
            //but note that the two indices must be successive, so maybe allocate
            //the nodes here and pass the index as parameter (c.f. child_index below)
            recurseBIH(b1, child_index, index_start, storeIndex, maxdepth-1);
            recurseBIH(b2, child_index+1, storeIndex, index_end, maxdepth-1);
        }
    }

    private bool intersectBIH_recurse(inout Ray isRay, inout Intersection inters,
        int nodeidx, inout Box b, bool visibility = false)
    {
        BIH_Node* node = &mBihNodes[nodeidx];

        uint axis = node.index_flags & BIH_Node.AXIS_MASK;
        //writefln("node= %d", node);

        version (OnDemandBIH) {
            //on demand calculation of BIH tree
            if (axis == 3 && node.count > cMaxTrisPerNode) {
                //writefln("huh! %d %d", nodeidx, node.count);
                recurseBIH(b, nodeidx, node.item, node.item+node.count, 1);
                axis = node.index_flags & BIH_Node.AXIS_MASK;
            }
        }

        uint next = node.index_flags >>> BIH_Node.INDEX_SHIFT;

        if (axis == 3) {
            return intersectNaive(isRay, inters, node.item, node.count, visibility);
        } else {
    version (NewBIH) {
            float tmin, tmax = void;
            //re-intersect
            int s1 = isRay.sign[axis];
            int s0 = 1-s1;
            Box b1 = b, b2 = b;
            b1.p[s0][axis] = node.clip[s1];
            b2.p[s1][axis] = node.clip[s0];

            Intersection inters2 = void;
            bool is_b1, is_b2;

            is_b1 = b1.intersectRay(isRay, tmin, tmax);
            if (is_b1) {
                is_b1 = intersectBIH_recurse(isRay, inters, next+s1, b1, visibility);
            }
            if (is_b1 && visibility)
                return true;
            if (node.clip[0]>node.clip[1] || !is_b1) {
                if (b2.intersectRay(isRay, tmin, tmax))
                    is_b2 = intersectBIH_recurse(isRay, inters2, next+s0, b2, visibility);
                else
                    is_b2 = false;
            }

            if (is_b1 && !is_b2) {
                return true;
            } else if (!is_b1 && is_b2) {
                inters = inters2;
                return true;
            } else if (is_b1 && is_b2) {
                if (inters.t > inters2.t)
                    inters = inters2;
                return true;
            }
    }
    version (OldBIH) {
            float tmin, tmax = void;
            //re-intersect
            bool res = b.intersectRay(isRay, tmin, tmax);

            float ca = isRay.start[axis] + isRay.dir[axis] * tmin;
            float cb = isRay.start[axis] + isRay.dir[axis] * tmax;

            Box b1 = b, b2 = b;
            b1.p[1][axis] = node.clip[0];
            b2.p[0][axis] = node.clip[1];

            /*inters.triIdx = cast(int)((isRay.dir*tmin).length*20);
            inters.t = (isRay.dir*tmin).length;
            return true;*/
            bool b0_ok = false;
            bool b1_ok = false;

            if (ca <= node.clip[0]) {
                //box0
                if (intersectBIH_recurse(isRay, inters, next, b1, visibility))
                    return true;
                b0_ok = true;
            }
            if (ca >= node.clip[1]) {
                //box1
                if (intersectBIH_recurse(isRay, inters, next+1, b2, visibility))
                    return true;
                b1_ok = true;
            }
            if (!b0_ok && cb <= node.clip[0]) {
                //box0
                if (intersectBIH_recurse(isRay, inters, next, b1, visibility))
                    return true;
            }
            if (!b1_ok && cb >= node.clip[1]) {
                //box1
                if (intersectBIH_recurse(isRay, inters, next+1, b2, visibility))
                    return true;
            }
    }
        }
        return false;
    }

    private bool intersectNaive(inout Ray isRay, inout Intersection inters,
        uint start, uint count, bool visibility = false)
    {
        Intersection iTmp = void;
        bool hitFlag = false;
        inters.t = isRay.tMax;
        for (int i = start; i < start+count; i++) {
            bool hit;
            version (IntersectWald) {
                hit = mTriAccel[mTriIndex[i]].intersectRay(isRay, iTmp);
            } else {
                hit = mTris[mTriIndex[i]].intersectRay(isRay, iTmp);
            }
            if (hit && iTmp.t >= 0 && iTmp.t <= isRay.tMax) {
                if (visibility) {
                    //inters.triIdx = 0;
                    return true;
                }
                if (iTmp.t < inters.t) {
                    inters = iTmp;
                    inters.triIdx = mTriIndex[i];
                }
                hitFlag = true;
            }
        }
        return hitFlag;
    }

    private bool intersectBIH(inout Ray isRay, inout Intersection inters,
        bool visibility = false)
    {
        return intersectBIH_recurse(isRay, inters, 0, mBoundingBox, visibility);
    }

    this() {

/*        mTris = new Triangle[0];
        mNorms = new Triangle[0];
        mHasVertNorm = new bool[0];
        mTriData = new void*[0];*/
    }

    //middle point of scene
    public Vector3f middle() {
        return (mBoundingBox.p[0] + mBoundingBox.p[1])/2.0f;
    }

    public float maxExtent() {
        return fmax(fmax(fabs(mBoundingBox.p[0].x - mBoundingBox.p[1].x),
            fabs(mBoundingBox.p[0].y - mBoundingBox.p[1].y)),
            fabs(mBoundingBox.p[0].z - mBoundingBox.p[1].z));
    }

    public Box boundingBox() {
        return mBoundingBox;
    }

    public void addObject(AbstractObject obj) {
        Triangle[] tris = obj.getTriangles();

        TriData[] tmp = new TriData[tris.length];
        foreach (int i, inout TriData dat; tmp) {
            dat.obj = obj;
            dat.orgIdx = i;
        }

        mTris ~= tris;
        mTriData ~= tmp;

        assert(mTriData.length == mTris.length);

        update();
    }

    public void beginUpdate() {
        mUpdating = true;
    }

    public void endUpdate() {
        mUpdating = false;
        update();
    }

    public void update() {
        if (!mUpdating) {
            //precalculate
            mTriAccel.length = mTris.length;
            for (int i = 0; i < mTriAccel.length; i++) {
                mTriAccel[i].calculate(mTris[i]);
            }

            //reset the BIH related stuff (let the permutation array initialize)
            mTriIndex.length = 0;

            calculateBoundingBox();
            buildBIH();
        }
    }

    public void clear() {
        mTris.length = 0;
        mTriAccel.length = 0;
        mTriData.length = 0;
        mTriIndex.length = 0;
        update();
    }

    //calculate mBoundingBox and mAvgObjSize
    private void calculateBoundingBox() {
        mAvgObjSize = Vector3f(0,0,0);
        if (mTris.length > 0) {
            mBoundingBox = Box(mTris[0].v[0],mTris[0].v[0]);
            for (int i = 0; i < mTris.length; i++) {
                Box objectBox = Box(mTris[i].v[0], mTris[i].v[0]);
                for (int j = 0; j < 3; j++) {
                    if (mTris[i].v[j].x < mBoundingBox.p[0].x) mBoundingBox.p[0].x = mTris[i].v[j].x;
                    if (mTris[i].v[j].y < mBoundingBox.p[0].y) mBoundingBox.p[0].y = mTris[i].v[j].y;
                    if (mTris[i].v[j].z < mBoundingBox.p[0].z) mBoundingBox.p[0].z = mTris[i].v[j].z;
                    if (mTris[i].v[j].x > mBoundingBox.p[1].x) mBoundingBox.p[1].x = mTris[i].v[j].x;
                    if (mTris[i].v[j].y > mBoundingBox.p[1].y) mBoundingBox.p[1].y = mTris[i].v[j].y;
                    if (mTris[i].v[j].z > mBoundingBox.p[1].z) mBoundingBox.p[1].z = mTris[i].v[j].z;

                    //awww duplicated code!!!11111111
                    if (mTris[i].v[j].x < objectBox.p[0].x) objectBox.p[0].x = mTris[i].v[j].x;
                    if (mTris[i].v[j].y < objectBox.p[0].y) objectBox.p[0].y = mTris[i].v[j].y;
                    if (mTris[i].v[j].z < objectBox.p[0].z) objectBox.p[0].z = mTris[i].v[j].z;
                    if (mTris[i].v[j].x > objectBox.p[1].x) objectBox.p[1].x = mTris[i].v[j].x;
                    if (mTris[i].v[j].y > objectBox.p[1].y) objectBox.p[1].y = mTris[i].v[j].y;
                    if (mTris[i].v[j].z > objectBox.p[1].z) objectBox.p[1].z = mTris[i].v[j].z;
                }
                mAvgObjSize += (objectBox.p[1] - objectBox.p[0]);
                mAvgObjSize = mAvgObjSize / 2.0f;
            }
        } else {
            mBoundingBox = Box(Vector3f(0,0,0),Vector3f(0,0,0));
        }
    }

    /// Intersect a ray with the scene
    /// Returns: true on hit, false otherwise
    /// Params:
    ///  inters = will be filled with details about intersection on hit
    ///           undefined if no hit
    public bool intersect(inout Ray isRay, inout RenderIntersection rinters,
        bool visibility = false)
    {
        float tmin, tmax;
        if (mBoundingBox.intersectRay(isRay,tmin, tmax)) {
            assert(tmin <= tmax);

            Intersection inters;
            bool hit = false;
            version(DisableBIH) {
                hit = intersectNaive(isRay, inters, 0, mTris.length, visibility);
            } else {
                hit = intersectBIH(isRay, inters, visibility);
            }

            if (hit) {
                rinters.inters = inters;
                rinters.hitpoint = isRay.start+isRay.dir*inters.t;

                //get the triangle
                rinters.objtriidx = mTriData[inters.triIdx].orgIdx;
                rinters.object = getObject(inters);

                //texture coordinates
                Vector2f[] coords = rinters.object.getTexCoords(rinters.objtriidx);
                rinters.texcoords = coords[0]+(coords[1]-coords[0])*inters.u
                    +(coords[2]-coords[0])*inters.v;

                TriangleNorms triNorms = rinters.object.getTriNormals(
                    rinters.objtriidx, mNormInterpolate);
                rinters.normal = triNorms.getNormal(rinters.inters.u,
                    rinters.inters.v);

                rinters.material = rinters.object.getTriMaterial(rinters.objtriidx);
                if (!rinters.material)
                    rinters.material = globalMaterial;

                rayDiff(isRay,rinters,coords,triNorms);

                return true;
            }
        }
        return false;
    }

    ///hit test function, returns true if ray hits anything within ray range
    public bool intersectShadowRay(Ray isRay) {
        float tmin, tmax;
        if (mBoundingBox.intersectRay(isRay,tmin, tmax)) {
            assert(tmin <= tmax);

            Intersection inters;
            version(DisableBIH) {
                return intersectNaive(isRay, inters, 0, mTris.length, true);
            } else {
                return intersectBIH(isRay, inters, true);
            }
        }
        return false;
    }

    private void calcLPlane(inout Vector3f out_n, inout float out_d,
        in Vector3f a, in Vector3f b, in Vector3f c, in Vector3f norm)
    {
        out_n = (norm | (c-b)).normal;
        //Paper: Coeffs normalized so that L*a = 1
        float s = 1.0f/(out_n*(a-b));
        out_n *= s;
        out_d = out_n * b;
    }

    private void rayDiff(Ray isRay, inout RenderIntersection rinters,
        Vector2f[] coords, TriangleNorms triNorms)
    {
        float t = rinters.inters.t;
        Vector3f n = rinters.normal;
        Triangle tri = mTris[rinters.inters.triIdx];

        //Update dP/dx for current intersection

        Vector3f pxtdx = isRay.rdif_dPx + t*isRay.rdif_dDx;
        Vector3f pytdy = isRay.rdif_dPy + t*isRay.rdif_dDy;

        float dtdx = (pxtdx*n) / (isRay.dir*n);
        float dtdy = (pytdy*n) / (isRay.dir*n);

        isRay.rdif_dPx = pxtdx + dtdx*isRay.dir;
        isRay.rdif_dPy = pytdy + dtdy*isRay.dir;

        //Calculate L-planes for barycentric normal/tex-coord interpolation

        Vector3f La_N, Lb_N, Lc_N;
        float La_d, Lb_d, Lc_d;
        calcLPlane(La_N, La_d, tri.v[0], tri.v[1], tri.v[2], n);
        calcLPlane(Lb_N, Lb_d, tri.v[1], tri.v[2], tri.v[0], n);
        calcLPlane(Lc_N, Lc_d, tri.v[2], tri.v[0], tri.v[1], n);

        float La_dpdx = La_N*isRay.rdif_dPx;
        float Lb_dpdx = Lb_N*isRay.rdif_dPx;
        float Lc_dpdx = Lc_N*isRay.rdif_dPx;

        float La_dpdy = La_N*isRay.rdif_dPy;
        float Lb_dpdy = Lb_N*isRay.rdif_dPy;
        float Lc_dpdy = Lc_N*isRay.rdif_dPy;

        //Calculate differential texture coordinates

        rinters.rdif_dTx = La_dpdx*coords[0] + Lb_dpdx*coords[1] + Lc_dpdx*coords[2];
        rinters.rdif_dTy = La_dpdy*coords[0] + Lb_dpdy*coords[1] + Lc_dpdy*coords[2];

        //Calculate differential normals

        Vector3f dndx = La_dpdx*triNorms.n[0] + Lb_dpdx*triNorms.n[1] + Lc_dpdx*triNorms.n[2];
        Vector3f dndy = La_dpdy*triNorms.n[0] + Lb_dpdy*triNorms.n[1] + Lc_dpdy*triNorms.n[2];

        float ndotn = rinters.normal*rinters.normal;
        rinters.rdif_dNx = (ndotn*dndx - (n*dndx)*n) / pow(ndotn,1.5f);
        rinters.rdif_dNy = (ndotn*dndy - (n*dndy)*n) / pow(ndotn,1.5f);
    }

    public AbstractObject getObject(inout Intersection inters) {
        return mTriData[inters.triIdx].obj;
    }

    public Vector3f getNormal(inout RenderIntersection rinters) {
        return rinters.object.getTriNormals(rinters.objtriidx,
            mNormInterpolate).getNormal(rinters.inters.u,rinters.inters.v);
    }

    public int triCount() {
        return mTris.length;
    }

    public bool normInterpolate() {
        return mNormInterpolate;
    }
    public void normInterpolate(bool val) {
        mNormInterpolate = val;
    }
}
